MICROCOPY  RESOLUTION  TEST  CHART 

national  Bureau  OE  StanoaRQS  -  |»63  -  A 


AD- A 145  516 


REPORT  NO.  NADC-84027-30 


AUTO-REGRESSIVE  MOVING  AVERAGE 
SPECTRAL  ESTIMATION  USING  MODIFIED 
AND  LEAST  SQUARES  MODIFIED 
YULE-WALKER  ESTIMATES 


Charles  D.  Pizzichello 

Sensors  and  Avionics  Technology  Directorate 
Naval  Air  Development  Center 
Warminster,  Pennsylvania  18974 


2  FEBRUARY  1984 


FINAL  REPORT 


Approved  for  Public  Release; 
Distribution  Unlimited 


DTIC 


ELECTE 
SEP  10  1984 


Prepared  For 

NAVAL  AIR  DEVELOPMENT  CENTER 
Warminster,  Pennsylvania  18974 


NOTICES 


REPORT  NUMBERING  SYSTEM  —  The  numbering  of  technical  project  reports  issued  by  the  Naval 
Air  Development  Center  is  arranged  for  specific  identification  purposes.  Each  number  consists  of 
the  Center  acronym,  the  calendar  year  in  which  the  number  was  assigned,  the  sequence  number  of 
the  report  within  the  specific  calendar  year,  and  the  official  2-digit  correspondence  code  of  the 
Command  Office  or  the  Functional  Directorate  responsible  for  the  report.  For  example:  Report 
No.  N ADC-7801 5-20  indicates  the  fifteenth  Center  report  for  the  year  1978,  and  prepared  by  the 
Systems  Directorate.  The  numerical  codes  are  as  follows: 


CODE 

OFFICE  OR  DIRECTORATE 

00 

Commander,  Naval  Air  Development  Center 

01 

Technical  Director,  Naval  Air  Development  Center 

02 

Comptroller 

10 

Directorate  Command  Projects 

20 

Systems  Directorate 

30 

Sensors  &  Avionics  Technology  Directorate 

40 

Communication  &  Navigation  Technology  Directorate 

50 

Software  Computer  Directorate 

60 

Aircraft  &  Crew  Systems  Technology  Directorate 

70 

Planning  Assessment  Resources 

80 

Engineering  Support  Group 

PRODUCT  ENDORSEMENT  —  The  discussion  or  instructions  concerning  commercial  products 
herein  do  not  constitute  an  endorsement  by  the  Government  nor  do  they  convey  or  imply  the 
license  or  right  to  use  such  products. 


APPROVED  BY: 


DATE:  /  /9t/C 


UNCLASSIFIED 


SECURITY  CLASSIFICATION  OF  THIS  FACE  fWhwi  D at*  En'ar'dj 


REPORT  DOCUMENTATION  PAGE 


I.  REPORT  NUMBER 

NADC-84027-30 


4.  title  (■««»  Submit) 

Auto-Regressive  Moving  Average 
Spectral  Estimation  Using  Modified  and 
Least  Squares  Modified  Yule-Walker  Estimates 


KE.AU  INSTRUCTIONS 
BEFORE  COMPLETING  FORM 


1.  RECIPIENT’S  CATALOG  NUMBER 


S.  TYPE  OF  REPORT  R  PERIOD  COVERED 


FINAL 


t.  performing  oro.  report  number 


7.  AUTHORS  •; 


Charles  D.  Pizzichello 


s-  performing  organization  name  ano  address 
Sensors  &  Avionics  Technology 
Directorate,  Naval  Air  Development  Center 
Warminster,  PA  18974 


II.  CONTROLLING  OFFICE  NAME  ANO  ADDRESS 

Naval  Air  Development  Center 
Warminster,  PA  18974 


<2.  REPORT  DATE 

2  February  1984 


IS-  NUMBER  OF  PAGES 

84 


4.  MONITORING  AGENCY  NAME  4  AOORCSSff/  dlllmrmnt  trom  Controlling  Olllem)  IS.  SECURITY  CLASS.  (ol  thio  report) 

Unclassified 


IS.  DISTRIBUTION  STATEMENT  (ol  Chi,  Report; 


Approved  for  public  release;  distribution  unlimited 


17.  DISTRIBUTION  STATEMENT  (ol  Ih #  okotracl  ontorod  In  Black  20.  II  dllloront  from  Report; 


IB.  SUPPLEMENTARY  NOTES 


Appendices  A,  B 


19.  KEY  WOROS  (CMtfmit  on  ravertt  old*  It  noemmoory  and  Identity  by  Mock  i 


Auto-Regressive  Moving  Average  (ARMA) 
Parametric  Spectral  Estimation 
Yule-Walker  Equations 


20.  ABSTRACT  fConffmi*  on  rovoro 


•ary  m d  Identify  by  Motk  mmkot) 


Many  popular  contemporary  spectral  estimation  methods  invoke  the  modeling  of  the  observa¬ 
tion  time  series  by  a  rational  transfer  function.  These  techniques  employ  the  approximation 
of  the  second-order  statistical  relationships  using  the  commonly  known  Yule-Walker  equations. 
Solving  these  equations,  one  obtains  the  parameter  estimates  of  the  hypothesized  rational  trans¬ 
fer  function  model.  These  parameter  estimates  represent  a  set  of  autocorrelation  time  lags  of  the 
observation  time  series  and  are  ideally  selected  to  optimize  the  model  being  considered.  Because 
of  the  optimization  criteria,  there  are  numerous  methods/techniques  which  estimate  parameters 


DO  i  jan  *71  1473 


COITION  OF  I  NOV  ••  IS  OBSOLETE 
5/N  0102-  LR- 014-4401 


UNCLASSIFIED 

SECURITY  CLASSIFICATION  OF  THIS  PAGE  I 


UNCLASSIFIED 

SECURITY  CLASSIFICATION  OF  THIS  RACE  fWfcm  Dma  K*tm*0 


20.  ABSTRACT  (Continued) 

that  have  been  proposed  in  the  literature.  The  objective  of  this  paper  is  to  investigate  two  of  these 
methods  and  report  their  respective  performance. 

The  most  general  rational  transfer  function  time  series  model  is  the  Auto-Regressive  Moving 
Average  (ARMA)  model.  This  model  consists  of  a  rational  transfer  function  containing  both  poles 
and  zeroes.  Of  the  numerous  techniques  available  in  estimating  the  ARMA  parameters,  there  are 
two  which  utilize  the  Modified  Yule-Walker  equations.  This  paper  investigates  both  the  Modified 
Yule-Walker  (MYW)  and  Least  Squares  Modified  Yule-Walker  (LSMYW)  methods  in  estimating 
the  ARMA  process.  The  MYW  method  estimates  the  ARMA  parameters  using  a  minimal  set  of 
Yule-Walker  equations.  In  contrast,  the  LSMYW  method  utilizes  the  parametric  estimation  of  an 
overdetermined  set  of  Yule-Walker  equations. 


NADC-84027-30 


PREFACE 

This  study  was  also  submitted  as  a  computer  simulation  project  in  partial  fulfillment  for  the 
requirements  of  a  graduate  course  while  working  towards  the  degree  Doctor  of  Philosophy  in 
Electrical  Engineering  at  the  University  of  Rhode  Island.  This  work  was  done  while  the  author 
was  supported  under  the  Advanced  Graduate  Study  Award  Program  sponsored  by  the  Naval  Air 
Development  Center,  Warminster,  Pennsylvania. 

The  author  acknowledges  the  contributions  of  those  individuals  who  supported  the  comple¬ 
tion  of  this  research.  Specifically,  he  acknowledges  the  assistance  of  Professors  Donald  W.  Tufts 
and  Steven  M.  Kay,  both  members  of  the  Department  of  Electrical  Engineering  at  the  University 
of  Rhode  Island,  Kingston,  Rhode  Island. 


The  author  of  this  report  is  located  at  the 
Naval  Air  Development  Center 
Warminster,  Pennsylvania  18974 


NADC-84027-30 


\  ABSTRACT 

'  \ 

Many  popular  contemporary  spectral  estimation  methods  invoke  the  modeling  of  the  observa¬ 
tion  time  series  by  a  rational  transfer  function.  These  techniques  employ  the  approximation  of  the 
second-order  statistical  relationships  using  the  commonly  known  Yule-Walker  equations.  Solving 
these  equations,  one  obtains  the  parameter  estimates  of  the  hypothesized  rational  transfer  function 
model.  These  parameter  estimates  represent  a  set  of  autocorrelation  time  lags  of  the  observation 
time  series  and  are  ideally  selected  to  optimize  the  model  being  considered.  Because  of  the  optimi¬ 
zation  criteria,  there  are  numerous  methods/techniques  which  estimate  parameters  that  have  been 
proposed  in  the  literature.  The  objective  of  this  paper  is  to  investigate  two  of  these  methods  and 
report  their  respective  performance. 

The  most  general  rational  transfer  function  time  series  model  is  the  Auto-Regressive  Moving 
Average  (ARMA)  model.  This  model  consists  of  a  rational  transfer  function  containing  both  poles 
and  zeros.  Of  the  numerous  techniques  available  in  estimating  the  ARMA  parameters,  there  are  two 
which  utilize  the  Modified  Yule-Walker  equations.  This  paper  investigates  both  the  Modified  Yule- 
Walker  (MYW)  and  Least  Squares  Modified  Yule-Waiter  (LSMYW)  methods  in  estimating  an  ARMA 
process.  The  MYW  method  estimates  the  ARMA  parameters  using  a  minimal  set  of  Yule-Walker 
equations.  In  contrast,  the  LSMYW  method  utilizes  the  parametric  estimation  of  an  overdetermined 
set  of  Yule-Walker  equations. 
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INTRODUCTION 

Spectral  estimation  has  received  much  attention  in  a  variety  of  applications  such  as  radar, 
sonar,  underwater  acoustics,  seismology  and  speech  over  the  last  twenty  years.  The  estimation  of 
the  power  spectral  density,  or  simply  the  spectrum,  is  usually  based  on  the  use  of  the  Fast  Fourier 
Transform  (FFT).  Until  the  introduction  of  the  FFT  in  1965,  the  classical  method  of  estimating  the 
power  spectrum  was  the  Blackman  and  Tukey  approach.  It  should  be  noted,  however,  that  spectral 
analysis  can  be  classified  into  two  categories,  namely  the  non-parametric  methods  and  the  para¬ 
metric  methods. 

The  non-parametric  methods  do  not  model  the  spectrum  but  estimate  the  power  spectrum  by 
performing  operations  on  the  observed  data  set  directly.  These  methods  are  the  classical  or  tradi¬ 
tional  techniques  in  estimating  the  power  spectrum.  Examples  of  these  techniques  include  the 
Blackman-Tukey,  Periodogram,  Modified  Periodogram  and  the  Averaged  Periodogram  spectral 
estimators.  The  approach  of  operating  on  the  observed  data  set  directly  without  modeling  is  compu¬ 
tationally  efficient  and  usually  produces  reasonable  results.  However,  for  certain  signal  classes,  there 
are  inherent  performance  limitations  when  a  spectrum  is  estimated  by  use  of  traditional  methods. 
Two  prominent  limitations  are  frequency  resolution  and  spectral  "leakage"  caused  by  windowing 
the  observed  data  set  when  the  data  are  processed  directly.  These  two  limitations  are  highlighted 
when  a  traditional  method  is  used  to  analyze  short  data  records.  In  general,  since  most  data  records 
are  short  in  nature,  alternative  spectral  estimation  techniques  must  be  considered. 

In  contrast  to  the  non- parametric  techniques  discussed  above,  the  parametric  methods  assume 
some  particular  rational  system  model  that  underlies  the  production  of  the  observation  time  series. 
The  parametric  methods  estimate  the  parameters  of  an  observed  data  set  and  then,  by  substituting 
these  estimates  into  the  theoretical  power  spectral  density  implied  by  the  model,  yield  the  com¬ 
puted  power  spectra.  These  are  the  modem  methods  for  estimating  power  spectra  and  alternatives 
to  processing  data  sets  directly.  Several  examples  of  the  models  used  in  implementing  these  modern 
methods  are:  the  Auto-Regressive  Moving  Average  (ARMA)  model,  the  Auto-Regressive  (AR) 
model,  and  the  Moving  Average  (MA)  model.  Each  of  these  modern  parametric  estimation  methods 
assume  the  observed  time  series  is  generated  from  a  particular  model  excited  by  white  gaussian 
noise. 

Of  the  three  models  mentioned  above,  the  most  general  rational  system  model  is  the  ARMA 
process.  An  ARMA  model  possesses  a  rational  transfer  function  containing  both  poles  and  zeroes. 
The  more  specialized  models  like  the  MA  (all  zero)  or  the  AR  (all  pole)  are  derived  from  the  general 
ARMA  model.  In  general,  the  ARMA  model  seems  to  be  much  more  effective  than  the  MA  and  AR 
models  in  modeling  an  observation  time  series  because  it  contains  both  poles  and  zeroes.  However, 
in  practice,  most  practitioners  favor  utilizing  the  specialized  MA  or  AR  models,  because  the  ARMA 
model  is  computationally  more  complex,  proving  to  be  inefficient  in  use.  This  has  stimulated 
considerable  activity  aimed  toward  generating  computationally  efficient  ARMA  modeling  al¬ 
gorithms.  Because  of  this  recent  activity,  several  new  ARMA  modeling  algorithms  have  been  pro¬ 
posed  and  presented  in  the  literature. 

The  objective  of  this  paper  is  to  investigate  two  of  these  more  recent  techniques  and  evaluate 
their  respective  performance.  These  two  techniques  are  the  Modified  Yule-Walker  (MYW)  and  Least 
Squares  Modified  Yule-Walker  (LSMYW)  approaches  to  estimating  an  ARMA  model.  Like  any  of 
the  parametric  methods  discussed  above,  the  ARMA  model  employs  the  approximation  of  second- 
order  statistical  relationships  which  utilize  the  commonly  known  Yule-Walker  equations.  In  the  case 
of  an  ARMA  time  series,  since  both  poles  and  zeroes  are  contained  in  the  rational  transfer  function 
model,  a  special  set  of  Yule-Walker  equations  should  be  evaluated.  These  special  equations  are 
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sometimes  called  the  Modified  or  Extended  Yule-Walker  equations.  To  evaluate  the  MYW  and 
LSMYW  approaches  in  estimating  an  ARMA  model,  these  Modified  Yule-Walker  equations  will 
be  utilized. 

The  MYW  and  LSMYW  approaches  in  estimating  an  ARMA  mode!  and  sub-optimal  techniques 
in  the  sense  of  practical  on-line  or  real-time  processing.  These  methods  estimate  an  ARMA  process 
by  a  three-step  procedure.  This  procedure  consists  of  estimating  the  AR  and  MA  parameters  separ¬ 
ately  rather  than  jointly  as  required  for  the  optimal  real-time  parameter  estimation.  The  procedure 
first  estimates  the  AR  parameters,  then  inverse  filters  the  original  time  series  by  the  AR  parameters, 
and  finally  estimates  the  MA  parameters.  This  procedure  intuitively  decouples  the  AR  parameters 
and  the  MA  parameters,  making  possible  the  estimation  of  each  parameter  easily  but  contributes 
to  its  sub-optimal  nature  as  a  real-time  processor.  Even  though  these  techniques  are  sub-optimal, 
they  still  make  the  computational  load  for  estimating  ARMA  parameters  more  realizable.  It  should 
be  noted  in  this  research,  the  MA  parameters  will  be  estimated  using  Durbin's  [1]  algorithm. 

THEORY 

As  discussed  in  the  previous  section,  an  ARMA  process  consists  of  a  rational  model  transfer 
function  containing  poles  and  zeroes.  The  order  of  the  numerator  designates  the  number  of  zeroes 
q,  and  the  order  of  the  denominator  designates  the  number  of  poles  defined  by  p.  Since  an  ARMA 
(p,  q)  process  assumes  an  ARMA  transfer  function  excited  by  white  gaussian  noise,  the  observation 
time  series  x(n)  is  depicted  in  Figure  1  as  the  response  or  output  of  a  linear  time  invariant  system. 


s(n)  - 

White  Noise 


Bq(z) 

A^li) 


x(n) 


ARMA  (p,  q) 

Transfer  Function 

Figure  1 

The  output  x(n)  can  be  written  as  [3j 

p  q 

x(n)  =  -  y  ^  a(k)  x(n-k)  +  ^  ^  b(k)  n(n-k)  (1) 

k=1  k=0 

The  model's  AR  parameters  (the  a(k)'s)  and  the  MA  parameters  (the  b(k)'s)  can  be  obtained 
directly  by  a  set  of  exact  autocorrelation  lags  given  by 

°-xx(n)  =  E  (x(n+k)  x*  (k))  (2) 

where  E  (  )  denotes  the  expected  value  operation  and  x*  denotes  the  complex  conjugate. 

Since  the  auto-correlation  lags  are  complex,  it  can  be  assumed  that  the  negative  lags  will  be  com¬ 
plex  conjugate  symmetric  (i.e.,  Rxx(-n)  =  R*xx(n>).  It  should  be  noted  that  throughout  this 

paper  whenever  negative  lags  are  required  this  property  will  be  automatically  assumed. 
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Therefore,  by  multiplying  both  sides  of  equation  (1 )  by  x*  (n-L)  and  taking  the  expectation 
/'elds 


p  q 

Rxx»->  *  -  ^  a<k>Rxx'l-k>  +  2  b<k,Rnx'L-IO 

k=1  k=0 


(3) 


where  Rnx(L-k)  =  E  (n(n)x*(n-k)) 
since  Rnx(k)  =  0,  for  k  >  0,  then  equation  (3)  can  be  written  as  (3) 


RXX(L) 


q 

a(k)Rxx(L-k)  +  b{k)Rnx(L-k) 

k=0 


for  L  =  0,  1 ,  2 . q 


P 

-  ^  a(k)Rxx(L-k) 
k=1 


(4) 


(5a) 


(5b) 


for  L  =  q+1,  q+2 


Equations  (5)  represent  the  Yule-Walker  equations  for  an  ARMA  process.  Careful  examination 
of  the  Yule-Walker  equations  for  an  ARMA  process  indicates  that  the  higher  order  terms  do  not 
contain  terms  of  the  MA  parameters.  It  seems  that  for  L  greater  than  q,  the  MA  parameters  are 
decoupled  from  the  AR  parameters.  These  higher  order  equations  have  been  called  the  Modified 
or  Extended  Yule-Walker  equations.  Expanding  the  Modified  Yule-Walker  equations  given  by 
equation  (5b)  yields 


Rxx(q)  Rxx^-1^  •••  Rxx^-P+f) 

a1 

-Rxx(q+1) 

Rxx<0+1>  Rxx^ 

a2 

. 

-Rxx(q+2> 

Rxx(q+p-1)  Rxx(q+p-2)  ...  Rxx(q) 

3P 

-Rxx(q+p) 

(6) 


in  maxtrix  notation,  they  can  be  written  more  compactly  as 
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The  previous  section  described  the  two  methods  being  investigated  in  this  paper  namely,  the 
MYW  and  LSMYW  techniques.  Since  both  these  methods  compute  the  AR  parameters  first,  then 
compute  the  MA  parameters  second,  the  theory  for  computing  the  AR  parameters  for  these  two 
methods  is  presented. 

The  MYW  method  computes  the  AR  parameters  by  using  a  "minimal"  set  of  Yule-Walker 
equations.  It  utilizes  the  generation  of  the  Modified  Yule-Walker  equations  given  by  equation  (7). 
Using  this  approach,  the  matrix  R^x  of  equation  (7)  will  be  Toeplitz,  square,  non-symmetric  and 

not  assured  to  be  either  positive-definite  or  non-singular.  In  order  to  solve  for  the  AR  parameters, 
a  simple  matrix  inversion  of  R^x  yields 

a  =  _Rxx  Rxx  18) 


The  LSMYW  method  computes  the  AR  parameters  by  using  an  "overdetermined"  set  of 
Yule-Walker  equations.  This  method  was  studied  by  James  A.  Cadzow  in  1979  [2] .  Cadzow's 
method  hypothesizes  an  overdetermined  parametric  approach  which  reduces  the  undesired  hyper¬ 
sensitivity  caused  by  the  minimal  set  of  Yule-Walker  equations.  Using  this  approach,  however, 
causes  the  Rxx  data  matrix  from  equation  (7)  to  be  rectangular  or  non-square.  This  requires  var¬ 
ious  least-squares  techniques  to  be  implemented  which  compute  the  pseudoinverse  of  the  rec¬ 
tangular  data  matrix  R^x  and  allow  the  solving  of  equation  (8)  for  the  AR  parameters.  The  over¬ 
determined  input  data  matrix  using  the  Modified  Yule-Walker  equations  may  be  defined  as 


R(q)  R(q-1) 

...  R(q-p+D 

a1 

-R(q+1) 

• 

a2 

-R(q+2) 

R(M-1)  R(M-2) 

R(M-p) 

aP 

-R(M) 

where  M»q+p 


(9) 


in  matrix  notation,  they  can  be  written  more  compactly  as 


RM a  RM 


(10) 


To  solve  equation  (10),  since  R^j  is  non-square,  a  least  squares  procedure  should  be  utilized.  De¬ 
fining  equation  (10)  in  the  least  squares  sense  yields 


RM  RM  5 ■ -RS 


(ID 


where  H  is  the  complex  conjugate  transpose. 

Solving  for  the  AR  Parameters  using  equation  (11)  gives 

a  *  -<rm  rm)  1  RM  rm 


(12) 
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Now  that  the  computation  of  the  AR  parameters  for  both  the  MYW  and  LSMYW  methods 
has  been  defined,  our  attention  will  now  focus  on  computing  the  MA  parameters.  Recall  that 
in  both  methods,  the  next  step  requires  the  filtering  of  the  original  observation  time  series  of 
the  AR  components.  This  allows  the  computation  of  the  MA  parameters  to  be  calculated  separ¬ 
ately  and  efficiently.  Therefore,  utilizing  a  finite-impulse-response  (FIR)  filter  and  setting  its 
filter  coefficients  to  the  AR  parameters,  the  original  observation  time  series  can  be  filtered  yield¬ 
ing  a  pure  MA  process.  This  all-zero  filter  will  be  defined  as 

P 

Y(n)  =  a(k)  x(n-k)  (13) 

k=0 

where  x(n-k)  is  the  original  observation  time  series  a(k)  is  the  AR  parameters. 

The  third  and  final  step  is  now  to  compute  the  MA  parameters  using  some  appropriate  tech¬ 
niques.  In  this  research,  Durbin's  algorithm  was  utilized  which  computes  both  the  MA  parameters 
and  the  variance  ($2)  associated  with  that  process. 

Once  the  parameters  of  an  ARMA  (p,q)  process  are  computed,  including  the  variance  (s2) 
of  the  MA  process,  the  power  spectral  density  can  be  computed  by  [3] 

Px(f)  *  lH(exp(j27rf)  l2Pn  (f)  (14) 

q 

=  S2  1  +  ^  ]  b(k)exp(-j27rfk)  2 
k=1 

P 

1  +  }  ]  a(k)exp(-j2irfk)  2 

k=1 

Thus  far,  all  the  computations  of  the  ARMA  process  have  been  derived  knowing  the  exact 
auto-correlation  lags  of  the  process.  However,  when  a  practical  implementation  is  simulated,  an 
estimate  of  the  auto-correlations  must  be  computed.  In  this  research,  the  unbiased  auto-correlation 
estimate  function  was  used  in  estimating  the  auto-correlation  lag  time  series.  This  estimate  auto¬ 
correlation  function  is  defined  as 

N-1-k 

R(k)  =  ^  ^  x*  (t)  x(t+k)  (15) 

t=0 

The  hat  ( A)  over  the  R  in  equation  (15)  denotes  an  estimate. 
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SIMULATION  PROCEDURE 

The  investigation  of  restimating  an  ARMA  process  utilizing  the  MYW  and  LSMYW  methods 
was  performed  using  specially  developed  software  on  an  International  Business  Machines 
Corporation  (IBM)  370  mainframe  central  computer  system.  The  FFT  and  matrix  inversion  com¬ 
putations  were  accomplished  using  the  International  Mathematical  and  Statistical  Library  (IMSL). 
The  specialized  plotting  routines  were  developed  using  the  COMPLOT  software  algorithms.  Figure  2 
depicts  the  overall  software/system  flow  diagram. 

The  procedure  for  evaluating  each  respective  method  was  to  first  generate  an  ARMA  model 
time  series  given  a  set  of  ARMA  coefficients.  After  generating  the  ARMA  time  series,  the  ARMA 
coefficients  were  estimated  using  both  techniques,  then  compared  with  the  original  set  of  given 
ARMA  coefficients.  Table  1  depicts  the  given  initial  coefficient  test  conditions  for  an  ARMA  (4,  2) 
process  for  four  unique  test  cases.  A  somewhat  better  presentation  and  easier  interpretation  of  these 
initial  test  conditions  can  be  achieved  by  calculating  the  pole/zero  locations  for  each  unique  test 
case.  The  initial  pole/zero  test  conditions  were  computed  and  are  listed  in  Table  2.  Examining 
Table  2  shows  that  Case  #1  and  Case  #2  have  moderately  strong  poles,  while  Case  =f3  and  Case  =4 
exhibit  very  strong  poles.  The  radius  is  very  close  to  the  unit  circle.  In  contrast,  with  respect  to 
zero  locations,  Case  -2  and  Case  -=3  are  strong  while  Case  #1  and  Case  #4  are  weak.  An  overview 
of  Table  2  yields  the  features  of  each  unique  test  case  shown  in  Table  3. 

These  features  play  an  important  role  in  generating  the  ARMA  (4,  2)  time  series  and  also  af¬ 
fect  performance  results  as  discussed  later. 

The  first  step  of  the  simulation  procedure  consisted  of  the  generation  of  random  White 
Gaussian  Noise  (WGN)  which  was  used  to  generate  the  ARMA  (4,  2)  time  series.  In  generating  the 
ARMA  (4,  2)  times  series,  the  impulse  response  of  an  infinite  pulse  response  (MR)  filter  with  the 
ARMA  (4,  2)  coefficients  was  computed.  Then  the  impulse  response  and  the  WGN  were  convolved 
to  generate  the  ARMA  (4,2)  time  series.  It  should  be  noted  when  implementing  an  HR  filter, 
in  theory  its  length  is  infinite;  however,  in  order  to  implement  an  HR  filter  in  practice,  only  a  fin¬ 
ite  filter  length  should  first  be  computed  to  some  finite  decimal  point  accuracy.  The  finite  length 
of  the  HR  filter  computed  in  this  research  was  computed  by 

rn  =  0.00002  (16) 

where  r  =  radius  of  pole 

n  =  finite  filter  length 

The  finite  length  is  very  dependent  on  the  closeness  of  the  pole  to  the  unit  circle.  Therefore, 
in  computing  the  finite  filter  length,  the  strongest  pole  radius  for  each  case  was  used.  For  Cases 
#1  and  #2,  n  was  equal  to  40,  while  for  Cases  #3  and  #4,  n  was  set  to  500. 

Another  important  point  of  interest  is  when  the  convolution  process  occurs,  the  impulse  re¬ 
sponse  sequence  is  flipped  about  the  origin  and  then  slides  and  multiplies  the  WGN.  This  flipping 
operation  causes  the  index  of  the  finite  filter  response  to  go  from  positive  to  negative. 

This  causes  the  initial  values  of  the  convolution  process  to  be  non-stationary  and  contain 
transient  spikes.  To  generate  a  stationary  ARMA  (4,  2)  time  series,  these  initial  time  series  values 
should  be  filtered  out.  The  number  of  points  filtered  out  or  thrown  away  should  be  greater  than 
the  finite  filter  length  of  the  HR  filter  response.  For  added  safety,  assuming  worst  case  condi¬ 
tions  for  all  test  cases,  the  first  1000  time  series  values  were  discarded  to  ensure  stationarity  in  all 
simulation  runs. 
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Figure  2.  Software  Simulation  Flow  Diagram 
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Table  1 .  Initial  Coefficient  Test  Conditions 


A1 

A2 

A3 

A4 

81 

mm 

CASE  #1 

-1.352 

1.338 

-.662 

.240 

-.2 

.04 

CASE  #2 

-1 .352 

1.338 

-.662 

.240 

-.9 

.81 

CASE  #3 

-2.7604 

3.8085633 

-2.6535 

.9238 

-.9 

.81 

CASE  #4 

-2.7604 

3.8085633 

-2.6535 

.9238 

-.2 

.04 

Table  2.  Initial  Pole/Zero  Test  Conditions 


POLE  #1 

POLE  tt2 

ZEROrrl 

RADIUS 

FREQUENCY 

RADIUS 

FREQUENCY 

RADIUS 

FREQUENCY 

CASE  #1 

.700 

0.125 

.699 

0.208 

.200 

0.166 

CASE  #2 

.700 

0.125 

.699 

0.208 

.900 

0.166 

CASE  *3 

.979 

0.109 

.981 

0.140 

.900 

0.166 

CASE  *4 

.979 

0.109 

.981 

0.140 

.200 

0.166 

NOTE:  ALL  POLE/ZERO  VALUES  LISTED  ABOVE  ALSO  HAVE  SYMMETRIC  COMPLEX 
CONJUGATE  VALUES 


Table  3 


POLES 

ZERO 

Case  #1 

Moderately  Strong 

Weak 

Case  #2 

Moderately  Strong 

Strong 

Case  #3 

Strong 

Strong 

Case  #4 

Strong 

Weak 
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After  the  ARMA  time  series  had  been  generated,  the  next  step  was  to  estimate  the  auto¬ 
correlation  of  the  generated  time  series.  In  this  research,  the  unbiased  auto-correlation  estimate 
function  was  computed  using  equation  (15)  from  the  previous  section.  Following  the  estimation 
of  the  auto-correlation  function,  both  the  M  YW  and  LSMYW  equations  could  now  be  generated 
from  the  correlation  estimates.  In  the  case  of  the  MYW  method,  the  minimum  number  of  correla¬ 
tion  time  series  lags  required  to  form  these  equations  is  p+q  (the  order  of  both  the  AR  and  MA  para¬ 
meters).  During  simulation  runs,  25  time  series  lags  were  computed.  When  the  LSMYW  equations 
are  computed,  an  overdetermined  set  of  correlation  time  series  lags  are  required.  This  length  must 
be  much  greater  than  p+q  and  is  designated  Ml.  During  LSMYW  simulation  runs,  Ml  was  set  at 
50,  and  75  time  series  lags  were  computed.  For  all  simulation  runs,  the  length  of  the  input  data  time 
series  was  256  points. 

The  next  step  was  to  estimate  the  AR  parameters  from  the  newly  formed  MYW  and  LSMYW 
equations.  Both  methods  require  a  matrix  inversion  which  provides  the  AR  parameter  estimates. 

In  the  LSMYW  approach,  an  intermediate  step  had  to  be  performed  prior  to  computing  its 
matrix  inversion.  This  additional  step  was  required  because  in  the  LSMYW  approach,  the  correla¬ 
tion  matrix  is  non-square,  therefore  requiring  a  pseudoinverse  of  the  correlation  data  matrix  which 
allowed  solving  the  AR  parameters. 

Next  the  original  time  series  was  filtered  by  a  FIR  filter  using  the  AR  parameters  as  filter  co¬ 
efficients.  This  filtering  operation  yielded  a  time  series  containing  only  MA  parameters.  Then  the 
MA  parameters  could  be  estimated  using  Durbin's  algorithm  which  includes  the  estimates  of  the 
MA  parameters  and  the  variance  associated  with  that  process.  The  final  step  was  to  compute  the 
power  spectral  density  using  the  estimated  AR  parameters,  MA  parameters  and  variance  and  plot 
the  results. 

RESULTS  AND  CONCLUSIONS 

The  results  of  the  simulations  of  the  MYW  and  LSMYW  methods  in  estimating  the  parameters 
of  an  ARMA  (4,  2)  process  are  contained  in  Appendix  A.  For  each  test  case  of  each  method,  there 
were  three  experiments  which  yielded  three  unique  plot  presentations.  These  three  presentations 
are  ( 1 )  the  true  power  spectrum  vs  the  estimated  spectrum  for  one  realization,  (2)  the  true  spectrum 
vs  an  average  of  50  estimated  realizations  and  (3)  the  multiple  plot  of  50  estimated  events.  In  addi¬ 
tion,  there  is  also  a  plot  of  the  true  spectrum,  a  true  pole/zero  plot  and  the  estimated  pole/zero  plot 
for  one  realization  for  each  test  case.  For  all  50  event  runs,  there  was  no  redundancy  or  overlap  of 
the  input  time  series  values  from  one  time-frame  to  the  next. 

Table  4  shows  the  estimated  coefficients  for  one  realization  of  the  LSMYW  method,  while 
Table  5  lists  the  estimated  coefficients  of  one  realization  of  the  MYW  method.  Comparing  the  re¬ 
sults  of  Table  4  and  5  with  Table  1,  indicates  that  in  general  the  LSMYW  method  is  superior  to  the 
MYW  method  in  estimating  the  parameters  of  an  ARMA  process. 

The  LSMYW  method  may  be  assumed  to  be  superior  to  the  MYW  method  because  the  LSMYW 
method  uses  an  overdetermined  set  of  correlation  lags.  For  the  cases  where  the  poles  are  close  to  the 
unit  circle,  the  LSMYW  method  should  outperform  the  MYW  method. 

A  careful  review  of  all  the  results  contained  in  Appendix  A  reveals  that  both  methods  exhibit 
hypersensitivity  caused  by  the  closeness  of  either  the  pole  or  zero  radius  to  the  unit  circle.  Cases  #3 
and  #4,  the  strong  pole  location  test  cases,  seem  to  exhibit  the  best  results  for  both  methods.  How¬ 
ever,  the  LSMYW  method  outperforms  the  MYW  method  as  expected.  In  Cases  #1  and  #2,  both 
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Table  4.  One  Realization  of  Least  Squares  Modified 
Yule-Walker  Estimates 


A1 

A2 

A3 

A4 

B1 

B2 

CASE  #1 

-1.8444 

2.1362 

-1.2857 

0.5029 

-0.6799 

0.2576 

| 

CASE  *2 

-1.2690 

1.0800 

-0.6429 

0.3318 

-0.7268 

0.3946 

CASE  #3 

-2.7511 

3.8050 

-2.6473 

0.9241 

-0.8457 

0.6776 

CASE  #4 

-2.7122 

3.7443 

-2.6036 

0.9207 

-0.1346 

0.0663 

Table  5.  One  Realization  of  Modified  Yule-Walker  Estimates 


A1 

A2 

— 

A3 

A4 

B1 

B2 

CASE  #1 

-0.1435 

1.4551 

-0.8115 

0.4746 

0.0340 

0.4338 

CASE  #2 

-1.2672 

1.4643 

-0.8598 

0.2371 

-0.6344 

0.8144 

CASE  #3 

-2.2047 

3.0046 

-2.1169 

0.9251 

-0.6242 

0.8997 

CASE  #4 

-2.7097 

3.5733 

-2.3706 

0.7619 

-0.2122 

-0.1136 
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methods  suffer  decreased  performance  which  is  attributed  to  the  weaker  pole  locations.  The  worst 
case  is  depicted  in  Case  #2  for  either  method  as  seen  in  the  multiple  50  realization  plot.  This  case 
is  the  moderately  strong  pole/strong  zero  test  case  which  exhibits  poor  variance  estimates  using 
Durbin's  algorithm.  This  poor  variance  estimate  is  attributed  to  the  strong  zero/moderately  strong 
pole  location  combination,  since  the  variance  estimates  for  the  other  test  cases  were  closer  to  the 
original  variance  used  to  generate  the  original  ARMA  (4,  2)  time  series. 

In  conclusion,  the  LSMYW  and  MYW  spectral  estimators  have  been  shown  to  produce  reli¬ 
able  ARMA  spectral  estimates  for  certain  signal  classes.  In  addition,  both  methods  contribute  to 
reducing  the  computation  burden  of  computing  an  ARMA  spectral  estimate  making  possible  its 
computation  more  realizable.  It  has  also  been  shown  that  the  LSMYW  method  outperforms  the 
MYW  method  but  still  suffers  the  problem  of  hypersensitivity  caused  by  the  close  proximity  of 
either  the  pole  or  zero  location  to  the  unit  circle.  Therefore,  it  is  recommended  that  future  re¬ 
search  address  the  technical  issue  of  hypersensitivity  and  also  investigate  another  technique  in  es¬ 
timating  the  processes'  variance  independent  of  Durbin's  algorithm. 
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APPENDIX  A 

PLOTTED  EXPERIMENTAL  RESULTS 
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NOTE: 


LEGEND  FOR  APPENDIX  A 


SYMBOL 

DESCRIPTION 

MYWE 

Modified  Yule-Walker  Equations 

LSMYWE 

Least  Squares  Modified  Yule-Walker  Equations 

dB 

decibels 

X 

poles 

0 

zeroes 

ON  ALL  GRAPHS  THAT  CONTAIN  BOTH  CONTINUOUS  AND 
DASHED  LINES,  THE  DASHED  LINE  REPRESENTS  THE  TRUE 
SPECTRUM 
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APPENDIX  B 

SAMPLE  COMPUTER  FROGRAMS 


B-1 


w 


on  non  o  on  non  noon 
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job  <n<f pizz •. riM£  =  <2. 30 - 

exec  F0RtVCG. REGION. S0»1300Kt LIB 1-IMSLD.L I B2-CPL0T 
//r0RT .SYSIN  DD  * 

C  THIS  PROGRAM  COMPUTES  THE  MTU  SPECTRAL  ESTIMATE 

DIMENSION  W <  20000 ) 

real*8  h(o::ooo) »A(o:30> .B<o:30) 

REAL *8  X (0520000) .Y< 0:20000) »X3< 0520000) 

REAL >8  X2<  0  5  2000)  »A2<0.‘  30)  .  Y3<  052000) 

C3HPLEX*16  R(0:500)fXl<0:fC  ) .AK100.20) .Rl( 100) 
COMPLEX* 1 i  WA< 10200) .UM 1C  .  - 
CCMPL£X*16  Cl ( 100 ) »C  <  20 »  20 ) 

REAL *8  XF  <  300 ) . AL (100.100) . POy < 100 ) . X7 ( 300 > 

REAL*8  AYU< 100.100) .R7< 101 > .PYU< 100) 

C*t ************************** ***********************<***•*<** 

READ  THE  INPUT  DATA  AND  PARAMETERS 
NUN-NUMBER  OF  RUNS  REQUESTED 
N-LENGTH  OF  THE  IMPULSE  RESPONSE  REQUESTED 
IP-ORDER  OF  THE  AUTO-REGRESS I UE  PARAMETERS 
IQ -ORDER  OF  THE  MOVING-AVERAGE  PARAHETERS 
A (  I  ) -THE  AUTO-REGRESS I UE  COEFFICIENTS 
6(1) -  THE  MOVING  AVERAGE  coefficients 

READ ( 1.4) .N.NUM 
WRITE <6. 4) .N.NUM 
4  format  f 13 » 12 ) 

' E A D <  1 .2  >  •  IP  .  IQ 
URITE(6»2> . IP. 13 

2  FORMAT (212) 

DO  13  K - 1 . IP 
READ < 1 .3  >  >A<K  > 

3  FORMAT (FI  1.7) 

13  CONTINUE 

DO  16  J-l.IQ 
READ (1.3) »  B ( J ) 

16  CONTINUE 


GENERATE  GAUSSIAN  WHITE  NOISE 

NPTS-16000 

XSEED-1 1 111 

VAR=1 . 

CALL  WGN  <  W . NPTS  » VAR . I  SEED ) 

STORE  THE  WHITE  NOISE  IN  X<I> 

DO  10  I-O.NPTS-1 
X(T)-W(I+1) 

10  CONTINUE 

ft********************************************************** 

B(0)-1 .0 

GENERATE  THE  IMPULSE  RESPONSE 
CALL  IMPULS(N.H.A.B.IP.IO) 

100  FORMAT  <4(3X.F9.4) ) 

GENERATE  ARMA  TIME  SERIES 

LX-NPTS 

LH=N 

LY-NPTS+N 

CALL  CONVOL <H.LH»X.LX»Y.LT) 

FILTER  THE  NON-STATIONARY  DATA  AND  SAVE  TIME  SERIES 
NS  •  LENGTH  OF  THE  NON-STATIONARY  VALUES  OF  H<N) 

NS  *  1000 
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DO  17  I=0.LY-NS-1 
X3< I )=Y< I+NS) 

17  CONTINUE 

X1<I>  =  ARMA  STATIONARY  TIHE  SERIES  (COMPLEX) 

X2( I )  =  ARMA  STATIONARY  TIME  SERIES  (REAL) 

IFL  =  NUMBER  OF  STATIONARY  POINTS  REQUESTED 
IFL  MUST  BE  .LT.  Y(I+NS) 

URITE (  2  » 200 >  » NUM 
NUM1=0 

DO  SOO  KK2= 1 » NUM 
IFL=256 

DO  20  1=0. IFL-1 
XI ( I )=X3( I+NUM1 ) 

X2( I )=X3< I+NUM1 ) 

20  CONTINUE 

COMPUTE  THE  AUTO-CORRELATION  FUNCTION 
SET  LR=LENGTH  OF  AUTO-CORRELATION  LAOS  REQUESTED 
FOR  MYUE  LR-IP+IQ 

LR=23 
LN=IFL 

CALL  AUTO(Xl.R.LR.LN) 

C  GENERATE  THE  MODIFIED  YULE-WALKER  EQUATIONS 

CALL  MODY (R.IP.IQ.A1.R1. IP) 

C  INVERT  THE  MYWE  AND  SOLVE 

CALL  LEQ2C< A1  .IP.100.R1.1.  100 .0 .WA.WK »  IER ) 

URITE ( 6.200) IER 
200  FORMAT (12) 

C  SAVE  THE  AR  PARAMETERS  ESTIMATES 

DO  43  1=1 .IP 
A2< I )=R1 ( I ) 

45  CONTINUE 

C  FILTER  THE  AR  PARAMETERS 

NPTS1=IFL 

WRITE <6.1 00 )»(A2(K)»K=1»IP) 

CALL  ARFIL<X2» IP.A2.Y3.LY3.NPTS1 ) 

C  COMPUTE  THE  MOVING  AVERAGE  PARAMETERS 

CALL  DURBIN( Y3 .L Y3 . IQ . SO . AL » POWOl ) 
URITE<6»100)»(AL(I»2).I=1»IQ) 

URITE<6. 150). POWOl 
130  FORMAT (Fll .7) 

URITE (2.700) .(A2(K).K=l»IP)»<AL(I»IQ).I“l.I0) .POWOl 
700  FORMAT (21F1 1.7) 

URITE<6.200) »KK2 
NUM1=NUM1+IFL 
800  CONTINUE 
STOP 
END 
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//PIZZA  JOB  (NYF101) » 'PIZZ' .TIME=<2»30> 

//  EXEC  F0RTVCG»REGI0N.G0=1500K»L'IB1=IMSLD»LIB2=CPLQT 
//FORT .SYSIN  DD  * 

C  THIS  PROGRAMS  COMPUTES  THE  LSMYU  SPECTRAL  ESTIMATE 

DIMENSION  W< 20000) 

REAL*8  H<0:2000> .AC0J30) .B<0:30> 

REALX8  X  <  0 : 20000 )  .  Y  <  0 : 20000  >  . X3 ( 0 : 20000 ) 

REALX8  X2<0t2000) »A2<0 530) ,Y3 <0)2000 ) 

COMPLEXXl 4  R < 0 : 500 > . XI <  0  J  500 ) . A1 < 1 00 . 20 ) . R 1 < 100  > 
C0MPLEXX14  UA< 10200). UK (100) 

C0MPLEXX14  C1(100).C<20.20) 

REALX8  XF  <  300 )  » AL< 100. 100) »POW< 100) .X7<300> 

REALX8  AYU< 100.100) .R7< 101 ) ,PYW< 100) 
CXXX*XXX*XXXXXXXXXXXXXXXXX*XX*X**X*XXXXXXXXXXXXXX ************ 
C  READ  THE  INPUT  DATA  AND  PARAMETERS 

C  NUM=NUMBER  OF  RUNS  REQUESTED 

C  N=LENGTH  OF  THE  IMPULSE  RESPONSE  REQUESTED 

C  IP=ORDER  OF  THE  AUTO-REGRESSIVE  PARAMETERS 

C  I Q= ORDER  OF  THE  MOVING-AVERAGE  PARAMETERS 

C  A<I)=rHE  AUTO-REGRESSIVE  COEFFICIENTS 

C  B( I ) =THE  MOVING  AVERAGE  COEFFICIENTS 

READ< 1.4) » N.NUM 
4  FORMAT (13.12) 

READ< 1.2). IP. IQ 

2  FORMAT  <  212 ) 

DO  15  K=1 . IP 
READ  <1.3)»A(K) 

3  FORMAT  <  FI 1 . 7 ) 

15  CONTINUE 

DO  14  J=1»IQ 
READ  <1.3).B(J) 

14  CONTINUE 

Cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

C  GENERATE  GAUSSIAN  WHITE  NOISE 

NPTS- 16000 
ISEED-11111 
VAR=1. 

CALL  WGN(W.NPTS.VAR.ISEED) 

C  STORE  THE  WHITE  NOISE  IN  X<I> 

DO  10  1=0. NPTS-1 
X<I)=W(I+1) 

10  CONTINUE 

cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

B<0)=1 .0 

C  GENERATE  THE  IMPULSE  RESPONSE 

CALL  IMPULS<N.H. A.B. IP. IQ) 

100  FORMAT  <4<3X»F9.4>> 

C  GENERATE  ARMA  TIME  SERIES 

LX=NPTS 
LH=N 

LY-NPTS+N 

CALL  CONVOL (H.LH.X.LX.Y.LY) 

C  FILTER  THE  NON-STATIONARY  DATA  AND  SAVE  TIME  SERIES 

C  NS  =  LENGTH  OF  THE  NON-STATIONARY  VALUES  OF  H<N> 

NS  =  1000 
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DO  17  I=0fLY-NS-1 
X3( I )*Y( I+NS ) 

17  CONTINUE 

XI (I)  =  ARMA  STATIONARY  TIME  SERIES  (COMPLEX) 

X2< I )  =  ARMA  STATIONARY  TIME  SERIES  (REAL) 

I PL  =  NUMBER  OF  STATIONARY  POINTS  REQUESTED 
IFL  MUST  BE  .LT.  (LY-NS-1) 

WRITE(2f200) fNUM 
NUM1=0 

DO  800  KK2=lfNUM 
IFL=236 

DO  20  I=Of IFL-1 
XI ( I ) =X3( I+NUM1 > 

X2( I )=X3( I+NUM1 ) 

20  CONTINUE 

COMPUTE  THE  AUTO-CORRELATION  FUNCTION 
SET  LR=LENGTH  OF  AUTO-CORRELATION  LAGS  REQUESTED 
FOR  LSMYUE  LR.GT.M1 

LR=75 
LN=IFL 
Ml=50 

CALL  AUTO(XIfRfLRfLN) 

C  GENERATE  THE  MODIFIED  YULE-UALKER  EQUATIONS 

CALL  M0DY(RfIP»IQ»A1»R1.M1) 

IR=M1— IQ 

CALL  INUR(A1fR1fIPfIRfC>C1) 

C  INVERT  THE  MYUE  AND  SOLVE 

CALL  LEQ2C(C>IP'100>C1'1»100'0>UA>UK>IER) 

WRITE ( 6  r  200  > IER 
200  FORMAT (12) 

C  SAVE  THE  AR  PARAMETERS  ESTIMATES 

DO  45  1=1 » IP 
A2< I )=C1 < I ) 

45  CONTINUE 

C  FILTER  THE  AR  PARAMETERS 

NPTS1=IFL 

WRITE (6f 100) » (A2(K) »K=1fIP) 

CALL  ARFIL( X2f IPfA2fY3fLY3fNPTS1) 

C  COMPUTE  THE  MOVING  AVERAGE  PARAMETERS 

CALL  DURBIN ( Y3  fLY3f IQfSOf ALfPOWOI ) 
WRITE(6f100)f(AL(If2)fI=1»IQ) 

WRITE (6f 150) fPOWOI 
150  FORMAT (F10 .7) 

WRITE(2f700)f(A2(K) fK=1fIP>  f ( AL( I f IQ) f 1  =  1 f IQ) fPOWOI 
700  FORMAT (2 IF 11. 7) 

WRITE(<Sf200)fKK2 
NUM1-NUM1+IFL 
800  CONTINUE 
STOP 
END 
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cxxxxxxxx ******************************************** ****** 

THIS  PROGRAM  GENERATES  WHITE  GAUSSIAN  RANDOM  NOISE 

INPUT  PARAMETERS 
WAR  =  THE  VARIANCE 
I SEED  =  11111 

NPTS  =  NUMBER  OF  DATA  POINTS  REQUESTED 

OUTPUT  PARAMETERS 
W(I>  =  THE  OUTPUT  TIME  SERIES 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

c 

SUBROUTINE  UGN<U»NPTSfVAR» ISEED) 

DIMENSION  W<512> 

M=NPTS 

IF  <MOD(NPTS»2) .NE.O)M=NPTS+l 
DO  10  1=1 *M 

CALL  RANDU ( ISEED* IY» YFL ) 

ISEEO=IY 
10  U ( I ) =YFL 

CALL  NORRAN<U*M.O. *VAR> 

RETURN 

END 

SUBROUTINE  NORR AN  <  W  *  M  * RME AN * VAR ) 

DIMENSION  U(512> 

PI =4 • OXATAN  < 1 • 0 ) 

L=M/2 

DO  10  1=1 fL 
U1=U  <  2X1-1 ) 

U2=W<  2X1 ) 

TEMP=SC1RT  ( -2  .  OXALOG  <  U1 )  ) 

U( 2X1-1 )  =  TEMPXCOS<  2 . 0XPIXU2) *SQRT ( VAR ) +RMEAN 
10  U  C  2X1 > =TEMPXSIN<  2 . 0XPIXU2 ) XSQRT ( VAR  >+RMEAN 

RETURN 
END 

SUBROUTINE  RANDU < IX . I Y * YFL ) 

I Y=IXXi553? 

IF< IY>5*4.6 

5  I Y= I Y +21 47483647+1 

6  YFL=IY 

YFL=YFLX . 4656613E-? 

RETURN 

END 
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c*******t**x***t********************x**»************«*(***t 

THIS  PROGRAM  CALCULATES  THE  IMPULSE  RESPONSE  OF  AN 
AUTO-REGRESS I UE/MOUING  AVERAGE  FILTER  TRANSFER  FUNCTION 

INPUT  PARAMETERS 

A ( I )  »  THE  AUTO-REGESSIVE  COEFFICIENTS 
B< I )  -  THE  MOVING  AVERAGE  COEFFICIENTS 
IP  =  THE  ORDER  OF  THE  AUTO-REGRESSIVE  COEFFICIENTS 
IQ  =  THE  ORDER  OF  THE  MOVING  AVERAGE  COEFFICIENTS 
N  »  NUMBER  OF  FILTER  RESPONSES  REQUESTED 

OUTPUT  PARAMETERS 
H  ( N )  =*  OUTPUT  DATA  VECTOR 

c»«***»****************s******x********y*******»**«******** 

c 

SUBROUTINE  IMPULS<N» H»A»B» IP» IQ) 

REAL *8  H ( 0 : 2000 )  *  A  ( 0 ! 30  ) » B  ( 0  S  30 ) 

H(0)«B<0) 

DO  10  I  =  1 » N—  1 
IC=I 

H( I ) =ODO 

IF  (I.GT.IP)  IC=IP 

DO  20  J=1»IC 

H( I ) =H( I ) -A  <  J ) <H< I— J ) 

20  CONTINUE 

IF  (I.GT.IQ)  GO  TO  10 
H( I ) *H< I ) +B < I  ) 

10  CONTINUE 
RETURN 
END 
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c********************************************************* 

C  THIS  PROGRAM  CONVOLVES  THE  INPUT  SIGNEL  X(N) 

C  WITH  THE  FILTER  IMPULSE  RESPONSE  H(N) 

C 

C  INPUT  PARAMETERS 

C  X(I)  =  THE  INPUT  TIME  SERIES 

C  LX  »  THE  LENGTH  OF  THE  INPUT  TIME  SERIES 

C  H<I>  =  THE  IMPULSE  RESPONSE 

C  LH  =  THE  LENGTH  OF  THE  IMPULSE  RESPONSE 

C 

C  OUTPUT  PARAMETERS 

C  r < I )  =  THE  OUTPUT  TIME  SERIES 

C  LY  =  THE  LENGTH  OF  THE  OUTPUT  TIME  SERIES 

C 

C  ****  LX  SHOULD  BE  GE.  LH  **** 

C 

C ******************************************* ************** 

c 

SUBROUT I NE  CONVOL  C  H  *  LH  »  X  »  LX  f  Y , L Y ) 

REAL*8  H(0 J2000) . X(0J2000) >Y (012000) 

INTEGER  LHfLX.LY 
DO  10  I=0»LY-1 
Y  <  I  )=>ODO 
N=MINO< IfLH-1) 

DO  20  J=O.N 
Y(I)»Y<I)+H(J)*X<I-J> 

20  CONTINUE 
10  CONTINUE 
RETURN 
END 
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CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

THIS  PROGRAM  CALCULATES  THE  UNBIASED  AUTO-CORRELATION 
ESTIMATE  OF  A  COMPLEX  FUNCTION 

INPUT  PARAMETERS 
LR  -  NUMBER  OF  LAGS  CALCULATED 
X1(I>  »  INPUT  TIME  SERIES 
LN  =  THE  LENGTH  OF  THE  INPUT  TIME  SERIES 

OUTPUT  PARAMETERS 

R< I )  a  THE  AUTO-CORRELATION  ESTIMATES 

cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

c 

SUBROUTINE  AUTO(Xl>R'LR»LN> 

INTEGER  LR 

C0MPLEXX16  R(0:300>»xi(0:s00) 

DO  10  K=0*LR-1 

N*LN-i-K 

R<K)"(ODO»ODO> 

DO  20  J=0»N 

R (K ) =R<  K ) +DCONJG ( XI <  J) ) XXI <  K+ J> 

20  CONTINUE 

M*LN-K 

R  <  K ) *R  <  K ) /FLOAT  <  M  > 

10  CONTINUE 

RETURN 
END 
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Cxyxxxxxxxyxxxyxxxxxyxyxxxxyxxxxxxxyyxxxxyxyxxxxxxyxyxxxxxyy 
THIS  PROGRAM  CALCULATES,  THE  MODIFED  YULE-UALKER 
ESTIMATES  BY  SORTING  THE  INPUT  AUTO-CORRELATION 
TIME  SERIES 

INPUT  PARAMETERS 

IP  *  THE  ORDER  OF  THE  AUTO-REGRESSIVE  PARAMETERS 
IQ  »  THE  ORDER  OF  THE  MOVING  AVERAGE  PARAMETERS 
R( I )  =  THE  INPUT  AUTO-CORRELATION  VALUES 

OUTPUT  PARAMETERS 

A1<I»J>  =  THE  LEFT  SIDE  MODIFIED  ESTIMATES 
R1 < I )  *  THE  RIGHT  SIDE  MODIFIED  ESTIMATES 

cxxxxxyyxyxxxxxxxxyxyyyyxxxxyxxyxyyyyxxxyxxxyxxxxxxyyxyxyxyx 

c 

SUBROUTINE  MODY<R > IP. 10. A1 »R1 .Ml > 

C0MPLEXX16  AK100.20)  .R<0!SOO)  .Rl(  tOO) 

IF ( IP.EQ.Ml )  M2=IP 
IF(IP.LT.Ml)  M2*Ml-IO 
IL*IP-1 
MaIG+1 

DO  10  1*1 .M2 
DO  20  J=0. IL 
KK=IQ— J+I— 1 

IF(KK.LT.O)  A1 < I . J+l ) *DCONJG<  R  < -KK ) ) 

IF(KK.GE.O)  A1(I.J+1)=R<KK) 

20  CONTINUE 
10  CONTINUE 

DO  30  1*1. M2 
K3=M+I-1 
R1 < I )=-R<K3) 

30  CONTINUE 
RETURN 
END 
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c»*  *************************************************  ******** 

THIS  SUBROUTINE  CONVERTS  A  NON-SQUARE  MATRIX  TO 
A  SQUARE  MATRIX  IN  ORDER  TO  SOLVE  ITS  INVERSE 

INPUT  PARAMETERS 

At ( I  *  J)  *  LEPT  SIDE  OP  NON-SQUARE  MATRIX 
R1(I)  =  RIGHT  SIDE  OP  NON-SQUARE  MATRIX 
IC  *  NUMBER  OP  COLUMNS  OF  THE  NON-SQUARE  MATRIX 
IR  -  NUMBER  OP  ROUS  OP  THE  NON-SQUARE  MATRIX 

OUTPUT  PARAMETERS 
C(I)  =  LEFT  SIDE  OP  SQUARE  MATRIX 
C1(I>  =  \IGHT  SIDE  OP  SQUARE  MATRIX 

C*********************************************************** 

c 

SUBROUTINE  INVRt A1 • R1 » IC » IR >C »C1 > 

C0MPLEXB16  Al<100>20)  *R1 < 100 > »C1 < 100>  »C(100»  100  > 

C  COMPUTE  A1  TRANSPOSE  R1 

DO  120  1=1 »IC 
Cl < I )=<ODOfODO> 

DO  120  J=1 t IR 

Cl ( I ) =C1 ( I )+DCONJG<  A1 < J> 1 ) ) *R1 <  J> 

120  CONTINUE 

C  COMPUTE  A1  TRANSPOSE  A1 

DO  130  1=1 fIC 
DO  130  J=1>IC 
C<I* J)=<ODO»ODO) 

DO  130  K=1 t IR 

C  < I » J) =C< I » J) +DCONJG ( A1 <K» I ) )*A1 <K » J> 

130  CONTINUE 
RETURN 
END 
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cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

FILTER  AR  PARAMETERS  AND  ESTIMATE  MA  PARAMATERS 

INPUT  PARAMETERS 
A2<  I )  =  AR  PARAMETERS 
IP  ■  ORDER  OF  AR  PARAMETERS 
X2<I>  *  ORIGINAL  INPUT  TIME  SERIES 
NPTS1  =  NUMBER  OF  INPUT  TIME  SERIES 

OUTPUT  PARAMETERS 
Y3<  I )  *  FILTERED  AR  TIME  SERIES 
LY3  =*  LENGTH  OF  FILTERED  TIME  SERIES 

cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

SUBROUTINE  ARFIL <X2» IP» A2» Y3 »LY3 rNPTSl > 

REALX8  A2<0130>  >X2<012000> .Y3<0£2000> 

M*NPTS1-IP-1 
LY3  =*  M+l 
A2C0>=1.0 
DO  10  NN»0*M 
Y3(NN+IP)=ODO 
DO  20  K-O.IP 

Y3<NN+IP)=Y3<NN+IP>+A2<K)XX2<NN+IP-K> 

20  CONTINUE 
10  CONTINUE 
RETURN 
END 
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c**************************************************** 

THIS  PROGRAM  CALULATES  THE  MA  PARAMETERS 

INPUT  PARAMETERS 

XF(N)=  INPUT  TIME  SERIES 

NPTSF*  NUMBER  OF  POINTS 

NQMAX*  ORDER  OF  THE  PROCESS 

NPL  *  NPTSF/3  (ORDER  FOR  LONG  AR  PROCESS) 

OUTPUT  PARAMETERS 

AL  (If  NOMAX  >  *  THE  B(I)  COEFFICIENTS 
POUOl  *  THE  VARIANCE 

**************************************************** 

SUBROUTINE  DURBIN <  XF f NPTSF » NQMAX . NPL r AL f POUOl ) 

REAL*8  XF  <  300) f AL( IOOf 100 ) fPOU< 100 ) 

CALL  LEVIN  <  XF f NPTSF f NPL fALf POUOl fPOU) 

XF(1>-1.0 
DO  10  1=1 » NPL 
10  XF(IM)*AL(IfNPL) 

NPL1=NPL+1 

CALL  LEV I N  <  XF . NPL 1 f NQMAX > AL » POUO » POU  > 

RETURN 

END 

SUBROUTINE  LEVIN( X7 . NPTS . NP  f AYUfPYUO t PYU ) 

THIS  SUBROUTINE  COMPUTES  THE  YULE-UALKER  ESTIMATES  OF  THE 
AR  PARAMETERS 

RE AL*8  X7 ( 300 ) . AYU <100»100)fR7<101)f PYU (100 ) 

N2=NP+1 
DO  20  K»1fN2 
NK=NPTS-K+1 
R7(K>=0.0 
DO  10  L*1 fNK 

R7<K)=R7(K)+X7(L)*X?(L+K-1) 

10  CONTINUE 

R7(K ) =R7 ( K ) /FLOAT ( NPTS ) 

20  CONTINUE 

AYU( 1 f 1 ) =-R7 ( 2 ) /R7 ( 1 ) 

PYU0=R7 ( 1 ) 

PYU(  1  )  =  (  1 .0-AYU(  1  f  1  )**2)*R7(  l ) 

IF ( NP  .EQ . 1 >G0  TO  60 
DO  50  I=2fNP 
11=1-1 
B*-R7< 1+1 > 

DO  30  K=1fI1 

30  B*B-AYU<Kf 1-1 ) *R7( I+l-K ) 

AYU(IfI)*B/PYU(I-1> 

DO  40  K*1 f II 

40  AYU(KfI)»AYU<KfI-1)+AYU(IfI)*AYU<I-KfI-1) 

PYU< I >=( 1 .0-AYU( IfI)**2) *PYU< 1-1 ) 

30  CONTINUE 

40  RETURN 

END 
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C********M*************************************************** 

THIS  SUBROUTINE  COMPUTES  THE  POWER  SPECTRAL  DENSITY 
OF  AN  ARMA  PROCESS  GIVEN  THE  MOVING  AVERAGE  AND 
AUTO-REGRESSIVE  COEFFICIENTS 

INPUT  PARAMETERS 

IQ  =  ORDER  OF  THE  MOVING  AVERAGE  PROCESS 
IP  *  ORDER  OF  THE  AUTO-REGRESSIVE  PROCESS 
A(I)  a  AUTO-REGRESSIVE  COEFFICIENTS 
B(I)  *  MOVING  AVERAGE  COEFFICIENTS 
N  *  LENGTH  OF  THE  FFT 
VAR  a  THE  VARIANCE  OF  THE  MA  PROCESS 

OUTPUT  PARAMETERS 
PS<I>  *  OUTPUT  POWER  SPECTRUM 

Ct*************t********************************************** 
SUBROUT INE  POWER  < 1 0  f I P  f  A  f  B  f  N  f  P  S  f  V AR  > 

R£AL*8  A(0:30) *B<0:30> fX5<2003 fX6<200) 

REAL<8  Y ( 200 ) • Y 1 < 200 )  f  PS  ( 200 ) 

C0MPLEX416  X8  <  200 ) f  X7 ( 200 ) 

REAL *8  WK(IOOO) 

DIMENSION  IWM1001) 

ND2*N/2 
N1*N— IP-1 
DO  10  I*OfIP 
X5< 1+1 ) *A< I ) 

10  CONTINUE 

DO  IS  1*1 fNI 
X5< 1+IP+l )*0D0 
13  CONTINUE 

CALL  FFTRC<X5fNfX7f IUKfWK) 

DO  20  I*2fND2 

X7 <  N+2- 1 ) -DCON JG  <  X7  < I >  > 

20  CONTINUE 

DO  23  1*1 fN 

X7  < I ) *DCON JG <  X7  < I ) > 

25  CONTINUE 

DO  26  1*1 fN 
Y<I)-CDABS<X7<1) >**2 

26  CONTINUE 

C  WRITE ( 6  f  60  >f(Y<K)fK*1fN> 

N2*N-IQ-l 
DO  30  I-0f10 
X6< I+l )*B< I ) 

30  CONTINUE 

DO  33  1*1 fN2 
X6(I+ia+l>*0D0 
33  CONTINUE 

CALL  FFTRC<X6fN»X8f IUKfWK) 

DO  41  I*2fND2 
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X8<  N+2-I )  =DCONJG (X3(  I )  ) 

41  CONTINUE 

DO  43  1*1  »N 

X8  < I > =DCON JG <  X8  < I >  > 

43  CONTINUE 

DO  44  1*1 »N 
ri<I)*CDABS<X8<I) )**2 

44  CONTINUE 

WRITE <6>60>. <Y1<K>  fK»lfN> 
DO  SO  1=1 fN 

PS<I)*<Y1<I)/Y<I) >*VAR 
SO  CONTINUE 

60  FORMAT  <5<3X»F?.6)> 

RETURN 

END 
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